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Abstract 

Boolean Delay Equations (BDEs) are semi-discrete dynamical mod- 
els with Boolean-valued variables that evolve in continuous time. Sys- 
tems of BDEs can be classified into conservative or dissipative, in a 
manner that parallels the classification of ordinary or partial differen- 
tial equations. Solutions to certain conservative BDEs exhibit growth 
of complexity in time. They represent therewith metaphors for bio- 
logical evolution or human history. Dissipative BDEs are structurally 
stable and exhibit multiple equilibria and limit cycles, as well as more 
complex, fractal solution sets, such as Devil's staircases and "fractal 
sunbursts." All known solutions of dissipative BDEs have stationary 
variance. BDE systems of this type, both free and forced, have been 
used as highly idealized models of climate change on interannual, in- 
terdecadal and paleoclimatic time scales. BDEs are also being used 
as fiexible, highly efficient models of colliding cascades of loading and 
failure in earthquake modeling and prediction, as well as in genetics. 
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In this paper we review the theory of systems of BDEs and illustrate 
their applications to climatic and solid-earth problems. The former 
have used small systems of BDEs, while the latter have used large 
hierarchical networks of BDEs. We moreover introduce BDEs with an 
infinite number of variables distributed in space ( "partial BDEs" ) and 
discuss connections with other types of discrete dynamical systems, in- 
cluding cellular automata and Boolean networks. This research- and- 
review paper concludes with a set of open questions. 

Keywords: Discrete dynamical systems, Earthquakes, El-Nino/Southern- 
Oscillation, Increasing complexity. Phase diagram. Prediction. 

PACS: 

02.30.Ks (Delay and functional equations) 

05.45.Df (Fractals) 

91.30.Dk (Seismicity) 

92. 10. am (El Nino Southern Oscillation) 

92.70.Pq (Earth system modeling). 

1 Introduction 

BDEs are a modeling framework especially tailored for the mathematical 
formulation of conceptual models of systems that exhibit threshold behavior, 
multiple feedbacks and distinct time delays [H [21 [3|, H] . BDEs are intended 
as a heuristic first step on the way to understanding problems too complex 
to model using systems of partial differential equations at the present time. 
One hopes, of course, to be able to eventually write down and solve the exact 
equations that govern the most intricate phenomena. Still, in the geosciences 
as well as in the life and other natural sciences, much of the preliminary 
discourse is often conceptual. 

BDEs offer a formal mathematical language that may help bridge the 
gap between qualitative and quantitative reasoning. Besides, they are fun 
to play with and produce beautiful fractals by simple, purely deterministic 
rules. Furthermore, they also provide an unconventional view on the concepts 
of non-linearity and complexity. 

In a hierarchical modeling framework, simple conceptual models are typi- 
cally used to present hypotheses and capture isolated mechanisms, while more 
detailed models try to simulate the phenomena more realistically and test for 
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the presence and effect of the suggested mechanisms by direct confrontation 
with observations [S]. BDE modehng may be the simplest representation 
of the relevant physical concepts. At the same time, new results obtained 
with a BDE model often capture phenomena not yet found by using con- 
ventional tools O d E]. BDEs suggest possible mechanisms that may be 
investigated using more complex models once their "blueprint" is detected 
in a simple conceptual model. As the study of complex systems garners in- 
creasing attention and is applied to diverse areas — from microbiology to 
the evolution of civilizations, passing through economics and physics — re- 
lated Boolean and other discrete models are being explored more and more 

[Sl[ini[IIl[I21[I3]. 

The purpose of this research-and-review paper is threefold: (i) summarize 
and illustrate key properties and applications of BDEs; (ii) introduce BDEs 
with an infinite number of variables; and (iii) explore more fully connections 
between BDEs and other types of discrete dynamical systems (dDS). There- 
fore, we first describe the general form and main properties of BDEs and 
place them in the more general context of dDS, including cellular automata 
and Boolean networks (Sect. [2]). Next, we summarize some applications, to 
climate dynamics (Sect. [3]) and to earthquake physics (Sect. Hj); these ap- 
plications illustrate both the beauty and usefulness of BDEs. In Sect. |5]we 
introduce BDEs with an infinite number of variables, distributed on a spatial 
lattice ("partial BDEs") and point to several ways of potentially enriching 
our knowledge of BDEs and extending their areas of application. Further 
discussion and open questions conclude the paper (Sect. [6]). 

2 Boolean Delay Equations (BDEs) 

BDEs may be classified as semi-discrete dynamical systems, where the vari- 
ables are discrete — typically Boolean, i.e. taking the values ("off") or 1 
( "on" ) only — while time is allowed to be continuous. As such they occupy 
the previously "missing corner" in the rhomboid of Fig. [1], where dynamical 
systems are classified according to whether their time [t) and state variables 
(x) are continuous or discrete. 

Systems in which both variables and time are continuous are called flows 
[T5] (upper corner in the rhomboid of Fig. [1]). Vector fields, ordinary 
and partial differential equations (ODEs and PDEs), functional and delay- 
differential equations (FDEs and DDEs) and stochastic differential equations 
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(SDEs) belong to this category. Systems with continuous variables and dis- 
crete time (middle left corner) are known as maps [TB| [T7] and include diffeo- 
morphisms, as well as ordinary and partial difference equations (OAEs and 
PAEs). 

In automata (lower corner) both the time and the variables are discrete; 
cellular automata (CAs) and all Turing machines (including real-world com- 
puters) are part of this group [TUl [HI [IB] , and so is the synchronous version 
of Boolean random networks fi2[ [T9] . BDEs and their predecessors, kinetic 
[20] and conservative logic, complete the rhomboid in the figure and occupy 
the remaining middle right corner. 

The connections between flows and maps are fairly well understood, as 
they both fall in the broader category of differentiable dynamical systems 
(DDS [m [151 HH])- Poincare maps ("P-maps" in Fig. 1), which are ob- 
tained from flows by intersection with a plane (or, more generally, with a 
codimension-1 hyperplane) are standard tools in the study of DDS, since 
they are simpler to investigate, analytically or numerically, than the flows 
from which they were obtained. Their usefulness arises, to a great extent, 
from the fact that — under suitable regularity assumptions — the process of 
suspension allows one to obtain the original flow from its P-map; hence the 
properties of the flow can be deduced from those of the map, and vice-versa. 

In Fig. 1, we have outlined by labeled arrows the processes that can lead 
from the dynamical systems in one corner of the rhomboid to the systems 
in each one of the adjacent corners. Neither the processes that connect the 
two dDS corners, automata and BDEs, nor these that connect either type of 
dDS with the adjacent-corner DDS — maps and flows, respectively — are as 
well understood as the (P-map, suspension) pair of antiparallel arrows that 
connects the two DDS corners. We return to the connection between BDEs 
and Boolean networks in Sect. 2.6 below. The key difference between kinetic 
logic and BDEs is summarized in Appendix A. 

2.1 General form of a BDE system 

Given a system with n continuous real- valued state variables v = {vi,V2, ■ ■ ■ , Vn) 
G for which natural thresholds G M exist, one can associate with each 
variable f j G M a Boolean-valued variable, G B = {0, 1}, i.e., a variable 
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I = 



n. 



(1) 



The equations that describe the evolution in time of the Boolean vector 
X = (xi,X2, . . . ,Xn) G B"" due to the time-delayed interactions between the 
Boolean variables Xj G B are of the form: 

= /i [t,Xi{t - 6'ii),X2(t - 612), ■ ■ ■,Xn{t - 9in)] , 
X2{t) = /2 [t,Xi{t - 92l),X2{t - 622), ...,Xnit~ 6'2„) , 

. (2) 

^ Xn(t) = fn [t,Xi{t - 9nl),X2(t - 9n2), ■ ■ ■ , Xn{t - 9nn)] ■ 

Here each Boolean variable Xi depends on time t and on the state of the 
other variables Xj in the past. The functions : B" — > B, 1 < i < n, are 
defined via Boolean equations that involve logical operators (see Table 1). 
Each delay value 9ij G M, 1 < ijj < n, is the length of time it takes for 
a change in variable xj to affect the variable Xj. One always can normalize 
delays % to be within the interval (0, 1] so the largest one has actually unit 
value; this normalization will always be assumed from now on. 

Following Dee and Ghil P, MuUhaupt [2], and Ghil and Mullhaupt, [3] 
we consider in this section only deterministic, autonomous systems with no 
explicit time dependence. Periodic forcing is introduced in Sect. [31 and ran- 
dom forcing in Sect. HI In Sects. 2-4 we consider only the case of n finite 
( "ordinary BDEs"), but in Sect. 5 we allow n to be infinite, with the variables 
distributed on a regular lattice ("partial BDEs"). 

2.2 Essential theoretical results on BDEs 

We summarize here the most important theoretical results from BDE theory; 
their original and complete form appears in [H [21 E]- 

We start by choosing a proper topology for the study of BDEs. Denoting 
by B"'[0, 1] the space of Boolean- valued vector functions with a finite number 
of jumps in the interval [0, 1]: 



and noting that, Vr, x |[t,t+i] still belongs to B"[0, 1] apart from a translation 
in time, the system can be considered as an endomorphism: 



x|[o,i]=x(t : 0<t< 1), 



J^f : B"[0,1] ^ B"[0,1]. 



(3) 
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We wish to extend this endomorphism into one that acts on the solutions 
x(t) of Eq. ([2]): 

■ X \ [t,t+l]^ X \ [t+l,t+2] ■ (4) 

Changing the point of view between ([2]) and (j4]) helps us study the dynamical 
properties of BDEs. The space B"[0, 1] equipped with Boolean algebra [21] 
and the topology induced by the Li metric: 

„1 n 

ci(x,y)= / y2\x,{t)-y,{t)\dt, (5) 

i=l 

is the phase space on which JF acts; we denote it by X. In coding theory, 
this metric is often called the Hamming distance. 

In constructing solutions for a given BDE system, there is a certain simi- 
larity with the theory of real-valued delay-differential equations (DDEs) (see 
[22l [23l [2 ^ [25] ) ■ as well as with that of ordinary difference equations (OAEs) 

mm)- 

Theorem 2.1 (Existence and uniqueness) Letx \[o,i]^ B"[0, 1] be the 
initial data of the dDS Q). Then the equivalent system ^ has a unique 
solution for all t > 1 and for an arbitrary n"^ -vector of delays 6 = {6ij) € 
(0, if. 

Sketch of Proof : The theorem can be proved by induction, constructing an 
algorithm that advances the solution in time and using a lemma that shows 
the number of jumps (between the values and 1) to be bounded from above 
in any finite time interval [1]. Thus the iterates JF'^, k = 1,...,K, stay 
within B"[0, 1] for all finite K and the unique solution of ([2]) is given simply 
by piecing together the successive intervals [0, 1], [1, 2], ... , [K, K + 1], etc. 
■ 

Theorem 2.2 (Continuity) The endomorphism T : X ^ X is con- 
tinuous for given delays. Moreover, the endomorphism JF : X x [0, 1]" 
X X [0, !]"■ is continuous, where the space of delays (0, 1]'" has the usual 
Euclidean topology. 

At this point, we need to make the critical distinction between rational 
and irrational delays. All BDE systems that possess only rational delays can 
be reduced in effect to finite cellular automata. Commensurability of the 
delays creates a partition of the time axis into segments over which state 
variables remain constant and whose length is an integer multiple of the 
delays' least common denominator (led). As there is only a finite number 
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of possible assignments of two values to these segments, repetition must 
occur, and the only asymptotic behavior possible is eventual constancy or 
periodicity in time. Thus, we obtain the following 

Theorem 2.3 ("Pigeon-hole" lemma) All solutions of ^ with ratio- 

2 

nal delays Q G Q" are eventually periodic. 

Remark. By "eventually" we mean that a finite-length transient may 
occur before periodicity sets in. An interesting feature of BDEs vs. flows or 
maps, as we shall see, is precisely that such transients have finite rather than 
infinite duration, i.e., asymptotic behavior is reached in finite time. 

Dee and Ghil pp, though, found that for the simple system of two BDEs: 



where v is the exclusive OR (see Table 1), the number of jumps per unit 
time seemed to keep increasing with time (see Fig. [2]) for a rational value 
6 = 0.977. Complex, aperiodic behavior only arises in cellular automata for 
an infinite number of variables (also called sites). Thus BDEs seem to pose 
interesting new problems, irreducible to cellular automata. One of these, at 
least, is the question of which BDEs, if any, do posses solutions of increasing 
complexity. To answer this question, we need to classify BDEs and to study 
separately the effects of rational and irrational delays. 



Based on the pigeon-hole lemma, and therefore on the behavior for rational 
delays, Ghil and Mullhaupt [3] classified BDE systems as follows. All sys- 
tems with solutions that are immediately periodic, for any initial data, are 
conservative] all other systems are dissipative and for some initial data will 
exhibit transient behavior before settling into eventual periodicity or quasi- 
periodicity. The DDS analogs are conservative {e.g., Hamiltonian) dynam- 
ical systems [28| [29] versus forced-dissipative systems {e.g., the well-known 
Lorenz system [30]). Typical examples of conservative systems occur in celes- 
tial mechanics [HT1I32], while dissipative systems are often used in modeling 
geophysical [331 [M] and many other natural phenomena. 

The simplest nontrivial examples of a conservative and a dissipative BDE 

are 



xi{t) 

X2{t) 



X2{t - 6) 

xi{t - 6) S7 X2{t - 1) 



(6) 



2.3 Classification of BDEs 



x{t) = x{t - 1) 
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and 

x{t) = x{t-i) Ax{t-e), o<e<i, 

respectively. The Boolean operators we use are listed in Table 1. It is com- 
mon to call a Boolean function / = (/i, . . . , /„) a connective and its argu- 
ments Xi channels [5T]; we shall also refer to a channel Xi simply as channel 
i. 

Definition 2.1 A BDE system is conservative for an open set Vt C 
(0, 1]'" of delays if for all rational delays in Vt and all initial data there 
are no transients; otherwise the system is dissipative. 

As is also the case in DDS theory, the conservative character of a BDE is 
tightly connected with its time reversibility. 

Definition 2.2 A BDE system is reversible if its time reversal also de- 
fines a system of BDEs. 

Theorem 2.4 (Conservative <^ Reversible) Definitions 2.1 and 2.2 
are equivalent. 

Useful algebraic criteria have been established |3] for linear or partially 
linear systems of BDEs to be conservative. Consider the following system 

n 

W = X] ^^^^ ~ ^'^^ ® 9i[xj'i't - I <i <n; (7) 

i=i 

here © and the summation symbol stand for addition (mod 2) in X; Cij G Z2, 
where Z2 is the field {0, 1} associated with this addition, while the Qi depend 
only on those Xj> for which Cij/ = 0. Note that x \j y = x (B y, while 
xAy = l©x©|/. We use the two types of symbols, y and ©, interchangeably, 
depending on the context or point of view. 

Adding constants Cjo to the above equations corresponds to adding par- 
ticular "inhomogeneous" solutions to the homogeneous linear system. All 
solutions of the full system can be represented as the sum of solutions to 
inhomogeneous and homogeneous systems. We review below only the homo- 
geneous case. 

We call a system linear ii and only if (iff) all Qi = 0. Naturally, the system 
obtained by putting (7^ = in ([7]) is called the linear part of the BDE system. 
Note that this concept of linearity (mod 2) is actually very nonlinear over 
the field of reals M, with usual addition and multiplication: it corresponds, 
in a sense, to the thresholding involved in Eq. ([1]). 

First we consider the simplest case of systems with distinct rational delays 
in their linear part. With any such system we associate its characteristic 
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polynomial 

Q{z) = det A{z), A,j = Sij + CijZ^'\ Vij = Q^tj, (8) 

where q is the led of all the delays 9ij such that Ctj ^ 0; the degree of Q is 
denoted by dQ. 

Theorem 2.5 (Conservativity for linear systems with distinct ra- 
tional delays) A linear system of BDEs is conservative for an open neigh- 
borhood Q of a fixed vector of distinct rational delays iff 

q sup 9ki = dQ. 

k 

In the case of rational delays only, we can give a first definition of partial 
linearity, namely that at least one Qi ^ and dQ > 2. 

Corollary 2.1 (Partially linear systems) The same result holds for a 
partially linear system of BDEs with distinct rational delays. 

2.4 Solutions with increasing complexity 

A natural question is whether (eventually) periodic solutions are generic in a 
BDE realm? We already noticed (see Fig. 2) that the answer to this question 
could be negative. Let us introduce the jump counting function 

J{t) = 7^{jumps of x(t) within the interval [t, t + 1)}, 

which measures the complexity of a BDE solution x(t) with a given set of 
initial data. 

Lemma 2.1 (Increasingly complex solutions for linear BDEs) All 

solutions ( except the trivial one x{t) = 0) of the linear scalar BDE 

x{t) = x{t - 1) V a;(t - ^2) V ■ ■ ■ V <t - O5) (9) 

with rationally independent < 9s < ■ ■ ■ < O2 < 9i = 1 and 6 > 2 are 
aperiodic and such that the lower bound for the corresponding J(t) increases 
with time. 

A simple example of this increasing complexity is given in Fig. [3|, for 
6 = 2,92 = 9 = (75-1)72, and a single jump in the initial data. Note that 
this delay is equal to the "golden ratio," which is the most irrational number 
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in the sense that its continued fraction expansion has the slowest possible 
convergence [35] , 

Remark. As for ODEs, a "higher-order" BDE can easily be written as 
a set of "first-order" BDEs ([2]). Therefore the previous lemma also applies 
to the system iQ of two linear BDEs, showing that the complexity of the 
solution is really increasing with time at least for irrational 6. 

A more general result holds for partially linear systems that include irra- 
tional delays. For such systems of the form ([7]), we introduce a generalized 
characteristic polynomial (GCP): 

g(A) = det A{\), Ay(A) = + CijX^'^. 

Clearly, this polynomial reduces to the characteristic polynomial in ([8]) if all 
the delays are rational and A = z'^. The index v of the GCP is defined as the 
number of its terms. We say that a BDE system ([7]) is partially linear if at 
least one (^j 7^ and v is large enough, z/ > 3. 

Theorem 2.6 (Increasingly complex solutions for partially linear 
BDEs) A partially linear system of BDEs has aperiodic solutions of increas- 
ing complexity, i.e. with increasing J{t), if its linear part contains 6 > 2 
rationally independent delays. 

The condition in this theorem is sufficient, but not necessary. A simple 
counterexample is given by the third-order scalar BDE 

x{t) = [x{t - 1) V xit - 6)] A x{t - r), (10) 

with 6, r, and 6/t irrational, and a single jump in the initial data at to- 
0<1 — ^<1 — r<to<l- The jump function for this solution grows in 
time like that of Eq. ([H]) for 5 = 2, although the GCP is identically 1, so that 
its index is = 1. 

On the other hand, there exist nonlinear BDE systems with arbitrarily 
many incommensurable delays that have only periodic solutions. For exam- 
ple, all solutions of 

x{t) = l[x{t-ek) (11) 

k=l 

are eventually periodic, with period tt = ^ 6'^ for n even, and vr = 2 ^ 6*^ 
for n odd; the length A of transients is bounded by A < vr. The multiplication 
in ffTTj) is in the sense of the field Z2, with xy = x Ay (see Table 1). 
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Dee and Ghil [T] established the upper bound on the jump function, 
J(t) < Kt''-\ where / is, in general, the number of distinct delays and the 
constant K depends only on the vector of delays 0. This bound is essential 
in proving the existence and uniqueness theorem in Sect. 12.21 Moreover, Ghil 
and MuUhaupt |3] obtained the lower bounds J{t) = O (ti°g2(^+i)) for Eq. 
and J{t) > K'£^°^'^^ for partially linear BDEs with 5 > v — 1 rationally 
independent delays in the linear part. These authors also showed the log- 
periodic character of the jump function in Fig. [3] (see also Fig. 7 in [3j). 

Having summarized these results, we are still left with the question why 
Fig. [2] here, with 9 = 0.997 being a rational number, does exhibit increasing 
complexity? The question is answered by the following "main approximation 
theorem" . 

Theorem 2.7 (Periodic approximation) All solutions to systems of 
BDEs can be approximated arbitrarily well (with respect to the Li-norm of 
X ), for a given finite time, by the periodic solutions of a nearby system that 
has rational delays only. 

The apparent paradox is thus solved by taking into account the length of 
the period obtained for a given conservative BDE and a given rational delay. 
As the led q becomes larger and larger, the solution in Fig. [3] here is well 
approximated for longer and longer times (see Fig. 9 of [3]); i.e., the jump 
function can grow for a longer time, before periodicity forces it to decrease 
and return to a very small number of jumps per unit time. 

Since the irrationals are metrically pervasive in M", i.e., they have mea- 
sure one, it follows that our chances of observing solutions of conservative 
BDEs with infinite — or, by the approximation theorem, arbitrarily long — 
period are excellent. In fact, the solution shown in Fig. 2 here was discovered 
pretty much by chance, as soon as Dee and Ghil [1] considered a conservative 
system. 

Ghil and Mullhaupt [3] studied, furthermore, the dependence of period 
length on the connective / and the delay vector G, as well as the degree of 
intermittency of self-similar solutions with growing complexity. In the latter 
case, we can consider each solution as a transient of infinite length. As we 
shall see next, such transients preclude structural stability. 
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2.5 Dissipative BDEs and structural stability 



The concept of structural stability for BDEs is patterned after that for DDS. 
Two systems on a topological space X are said to be topologically equivalent 
if there exists a homeomorphism h : X ^ X that maps solution orbits from 
one system to those of the other. The system is structurally stable if it is 
topologically equivalent to all systems in its neighborhood [T^ 

In discussing structural stability, we are interested in small deformations 
of a BDE leading to small deformations in its solution. A BDE can be 
changed by changing either its connective / or its delay vector O. Changes 
in / have to be measured in a discrete topology and cannot, therefore, be 
small. It suffices thus to consider small perturbations of the delays. 

Theorem 2.8 (Structural stability) A BDE system is structurally 
stable iff' all transients and all periods are bounded over some neighborhood 
U C M" of its delay vector 0. 

The periodic approximation theorem (Theorem 2.7) implies that, for 
BDEs like for DDSs, conservative systems are not structurally stable in 
X X [0, 1]" . Moreover, the conservative "vector fields," here as there, are in 
some sense "rare"; for BDEs they are just the three connectives x, x \/ y, 
and xAy, for which the number of O's equals the number of I's in the "truth 
table." Incidentally, the jump set on the delay lattice (see Figs. 1 and 3 
of [3]), and hence the growth of J{t), is exactly the same when replacing 
f{x,y) = xsyyhy f{x,y) = x A y |3]. 

The structural instability and the rarity of conservative BDEs justifies 
studying in greater depth dissipative BDEs. Ghil and MuUhaupt [3] concen- 
trated on the scalar rath-order BDE 



The connective / is most conveniently expressed in its normal forms from 
switching and automata theory, with xy = x A y and x + y = x V y. With 
this notation, the disjunctive and conjunctive normal forms represent / as a 
sum of products and a product of sums, respectively. This formalism helps 
prove that certain BDEs of the form fll2p lead to asymptotic simplification, 
i.e., after a finite transient, the solution of the full BDE satisfies a simpler 
BDE. An illustrative example is 



x{t)=f[x{t-9,),...,x{t-9n)]. 



(12) 



x{t) 



x{t-9i)x{t-92). 



(13) 
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where either 9i or 62 can be the larger of the two. Asymptotically, the 
solutions of Eq. (IT^ are given by those of a simpler equation 

x{t) = x{t - Oi). 

Comparison with the asymptotic behavior of forced-dissipative systems 
in the DDS framework shows two advantages of BDEs. First, the asymptotic 
behavior sets in after finite (rather than infinite) time. Second, the behavior 
on the "inertial manifold" or "global attractor" here can be described explic- 
itly by a simpler BDE, while this is rarely the case for a system of ODEs, 
FDEs, or PDEs. 

Finally, one can study asymptotic stability of solutions in the Li-metric 
of X. We conclude this theoretical section by recalling that, for < 6* < 1 
irrational, the solutions of 

x{t) = x{t - 9) x{t - I) 

are eventually equal to x(t) = 0, except for x{t) = 1, which is unstable. 
Likewise, for 

x{t) = x{t - 9) + x{t - 1), 

x{t) = 1 is asymptotically stable, while x{t) = is not. More generally, one 
has the following 

Theorem 2.9 Given rationally unrelated delays G = {9k}, the BDE 

n 

x{t) = Y[x{t-9k) 

k=l 

has x{t) = as an asymptotically stable solution, while for the BDE 

n 

x{t) = J2<t-^k), 

k=l 

x{t) = 1 is asymptotically stable. 

To complete the taxonomy of solutions, we also note the presence of 
quasi-periodic solutions; see discussion of Eq. (6.18) in Ghil and MuUhaupt 

Asymptotic behavior. In summary, the following types of asymptotic 
behavior were observed and analyzed in BDE systems: (a) fixed point — the 
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solution reaches one of a finite number of possible states and remains there; 
(b) limit cycle — the solution becomes periodic; (c) quasi-periodicity — 
the solution is a sum of several incommensurable "modes"; and (d) growing 
complexity — the solution's number of jumps per unit time increases with 
time. This number grows like a positive, but fractional power of time t [Il,[2], 
with superimposed log-periodic oscillations [5]. 

2.6 BDEs, cellular automata (CAs) and Boolean net- 
works 

We complete here the discussion of Fig. 1 about the place of BDEs in the 
broader context of dynamical systems in general. Specifically, we concentrate 
on the relationships between BDEs and other dDS, to wit cellular automata 
and Boolean networks. 

The formulation of BDEs was originally inspired by advances in theoret- 
ical biology, following Jacob and Monod's discovery [H7] of on-off interac- 
tions between genes, which had prompted the formulation of "kinetic logic" 
[20l [38l |39] and Boolean regulatory networks fl2[ [T9l H2] . In the following, 
we briefly review the latter and discuss their relations with systems of BDEs, 
whereas kinetic logic is touched upon in Appendix A. 

In order to understand the links between BDEs and Boolean regulatory 
networks it is important to start by recalling some well known deflnitions 
and results about cellular automata (CAs), which were introduced by von 
Neumann already in the late 1940s [18]. Doing so here will also facilitate the 
discussion of our preliminary results on "partial BDEs" in Sect. [51 

One deflnes a CA as a set of Boolean variables {xj : i = 1, . . . , N} 
on the sites of a regular lattice in D dimensions. The variables are usu- 
ally updated synchronously according to the same deterministic rule Xi(t) = 
f[xi{t — 1), ... , XN^t — 1)]; that is the value of each variable Xi at epoch t is 
determined by the values of this and possibly some other variables {xj} at 
the previous epoch t — 1. In the simplest case of -D = 1 {i.e., of a 1-D lattice) 
and flrst-neighbor interactions only, there are 2^ possible rules / : ^ B, 
which give 256 different elementary CAs (ECAs) studied in detail by Wolfram 
[TTl Bo] . For a given /, they evolve according to: 

Xi{t) = f[xi^,{t-l),x,{t-l),Xi+i{t-l)], l<t<N, (14) 

For a flnite size A^, Eq. iHM is a particular case of a BDE system ([2|) with 
connective fi = f for all i and a single delay 6ij = 1 for all i and j (see 
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also Sect. One generally speaks of asynchronous CAs when variables 
at different sites are updated at different discrete times according to some 
deterministic scheme. Such asynchronous CAs still belong to a restricted 
class of BDEs with integer delays 9ij G N. 

When both the space and the time are discrete, a finite-size CA will 
ultimately display either a fixed-point or periodic behavior. An important 
advantage of the great simplicity of EGAs is that it allows for systematic 
studies and helps understand their behavior in the limit of iV — >^ oo. It 
can be shown that different updating rules can lead to very different long- 
time dynamics. Wolfram [lOlBT] divided EGAs into four universality classes, 
according to the typical behavior observed for random initial states and large 
sizes N: For rules in the first class the system evolves towards a fixed point. 
For rules in the second class the dynamics can attain either a fixed point or 
a limit cycle, but in this case the period is usually small and it remains small 
for increasing iV-values. For rules in the third class, though, the period of the 
limit cycle usually increases with the size N and it can diverge in the limit 

^ oo, leading to "chaotic" behavior. Finally, GAs in the fourth class 
are capable of universal computation and are thus equivalent to a Turing 
machine. 

A first generalization of GAs are Boolean networks, in which the Boolean 
variables {xj : i = 1,2, . . . , N} are attached to the nodes (also called ver- 
tices) of a (possibly directed) graph and they evolve synchronously according 
to deterministic Boolean rules, which may vary from node to node. A fur- 
ther generalization is obtained by considering randomness, in the connections 
and/or in the choice of updating rules. In particular, the NK model intro- 
duced by Kauffman [191 112] , is among the most extensively analyzed random 
Boolean networks (RBNs). This model considers a system of Boolean vari- 
ables such that each value Xi depends on K randomly chosen other variables 
Xj through a Boolean function drawn randomly and independently from 2^ 
possible variants. The connections among the variables and the updating 
functions are fixed during a given system's evolution, and one looks for aver- 
age properties at long times. Since the variables are updated synchronously, 
at the same discrete t-values, the evolution will ultimately reach a fixed point 
or a limit cycle for any given configuration of links and rules. 

Kauffman [191 112] proposed such NK RBNs as models of a regulatory 
genetic network, with different nodes corresponding to different genes. The 
activity of a gene Xi is regulated by the activity of the other K genes to which 
Xi is connected. The different attractors, whether fixed point or limit cycle. 
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are related to different gene expression patterns. In this interpretation, a 
limit cycle, i.e. a recurrent pattern, corresponds to a cell type and the 
period is that of the cell cycle. 

The NK model was initially studied for a uniform distribution of the 
updating rules. In this situation, for small K values, one finds on average 
a small number of fixed points and limit cycles. The lengths of the possible 
attractors remain finite in the limit of ^ oo, and therefore the network 
dynamics appears "ordered." For large K values, though, the model displays 
"chaotic" behavior; in this case, the average number of attractors as well as 
their average length diverge with and the difference between two almost 
identical initial states can increase exponentially with time. Furthermore, one 
observes a transition in parameter space between typical dynamics that is 
characterized by a large connected cluster of frozen variables and the opposite 
one with small separated clusters of frozen variables. The critical values of the 
parameters corresponding to this passage from an "ordered" to a "chaotic" 
regime can be evaluated by looking at the evolution of the Hamming distance 
between two trajectories that start from slightly different configurations [121 
HU US]. In particular, for a uniform distribution of the Boolean updating 
functions, the NK model is "critical" when K = 2. 

Kauffman [19] suggested that natural organisms could lie on or near the 
borderline between these two different dynamical regimes, i.e. "at the edge 
of chaos," where the system is still sufficiently robust against small perturba- 
tions but at the same time close enough to the chaotic regime to feel the effect 
of selection. Accordingly, a lot of attention has been devoted to the study 
of such critical networks, which can also be obtained for K > 2 with appro- 
priate, and possibly more realistic, choices of the updating rule distribution. 
Similar edge-of-chaos suggestions have been made in other applications of 
dynamical systems theory, including DDS and celestial mechanics [IB]. 

For large A^-values, even the problem of determining the fixed points of 
a generic regular RBN, with Ki = K for all is highly nontrivial. In the 
context of the modeling of genetic interactions, the solution to this problem is 
thought to represent different accessible states of the cell, possibly triggered 
by external inputs [1?]. This problem has been recently reformulated [ITJIIS] 
in terms of the zero-energy configurations of an appropriate Hamiltonian. 
In this formulation, statistical mechanics tools from spin-glass physics can 
be brought to bear on the problem; these tools have also been successfully 
extended to general optimization issues |49j . 

Irregular RBNs, in which the number of inputs is also a node- dependent 
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random variable, are obviously harder to analyze. There is increasing evi- 
dence [20] that many networks arising in very different natural contexts are 
"scale free," i.e. their node- dependent connectivity Ki is distributed ac- 
cording to a power law P{Ki) oc Ki~^ . This seems to be true as well for 
the distribution of the input connections of some genetic networks. In the 
irregular case, too, one still observes "critical" dynamical behavior, given 
a suitable distribution of the updating Boolean functions. Kauffmann and 
colleagues [51] have recently studied the stability properties of regular and 
irregular RBNs and their dependence on the distribution of connectivity Ki 
and/or Boolean functions. Inverse problems, in which one tries to deter- 
mine the Boolean rules leading to a particular type of behavior, have been 
considered in [32] . 

The dependence of the average number m of attractors and of their period 
length on the size N in critical RBNs is still a matter of debate. This issue 
is particularly relevant for genetic-network modeling, since the behavior of 
rh{N) is expected to be related [T91B3] to the number of cell types which are 
present in an organism characterized by a given number N of genes. Recent 
results [S5| EH ES] show that rh{N) increases faster than any power of N in 
regular synchronous RBNs, in which all variables are updated in parallel at 
the same discrete epochs. These findings suggest that using different updat- 
ing schemes could lead to more realistic behavior, with a slower increase of 
fh{N). Moreover, assuming that all the variables act synchronously, i.e. that 
they "move in lock-step," may be too drastic a simplification for correctly 
modeling a number of natural systems and, in particular, interacting genes 
[56] . In order to overcome this simplification, asynchronous Boolean networks 
with different updating procedures have been studied recently [57]. A recent 
review [SS] underlines the links of such generalized RBNs with asynchronous 
CAs and with "kinetic logic." 

From the point of view of connections with BDEs, the updating scheme 
introduced by Klemm and Bornholdt [58] is of particular interest; they con- 
sider a "critical" regular RBN with K = 2 and weakly fluctuating delays in 
the response of each node. The number of stable attractors in this system 
increases more slowly with system size than for synchronous updating. It 
seems therefore that RBNs in continuous time may be more realistic and 
may exhibit new and possibly unexpected types of behavior. 

Oktem et al. [59] have recently applied a BDE approach to a Boolean net- 
work of genetic interactions with given architecture. In this case continuous 
time delays are introduced according to the BDE formalism of [H [2], El H] . As 
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a result, more complicated types of behavior than in synchronously updated 
Boolean networks have been observed, and the dynamics of the system seems 
to be characterized by aperiodic attractors. 

In both |58] and [59] , the authors introduced a minimal time interval be- 
low which changes in a given variable are not permitted. Such a cutoff, or 
"refractory period" [22], may have a physiological basis in genetic applica- 
tions, but it rules out the presence of solutions with increasing complexity. 
Therefore, for a finite number of variables, this restriction must result in an 
ultimately periodic behavior; the asymptotic period, though, could be much 
larger than the one obtainable with usual Boolean networks, especially when 
considering conservative connectives and irrational delays. From this point 
of view, the implementation of continuous time delays in [SHI EH] is differ- 
ent than in BDEs [I]-[l] and is similar to the one adopted in "kinetic logic" 
[20l [38| [39] , whose precise connections with our formalism are discussed in 
Appendix A. 

In the applications of BDEs that we review in the next sections, one finds 
different mechanisms that lead to aperiodic solutions of bounded complexity, 
without the need of a cutoff; one could thus explore the possibility of similar 
behavior in genetic-interaction models as well. Summarizing, one can say 
that kinetic logic and the recently proposed genetic network models [58| [59]. 
as well as others recent generalizations of RBNs with deterministic updating 
|56j . can be viewed either as asynchronous CAs or as particular cases of BDEs 
with large A^. 

In Sect. 5, we initiate the systematic study of BDE systems in the limit of 
an infinite number of variables, assumed for the moment to lie on a regular 
lattice and to interact according to a given, unique, deterministic rule. This 
study should allow us to better understand the connections of BDEs with 
(infinite) CAs, on the one hand, and with PDEs on the other. Such a study 
should also help clarify further the behavior of, possibly random. Boolean 
networks in continuous time. 

We now turn to an illustration of BDE modeling in action, first with 
a climatic example and then with one from lithospheric dynamics. Both of 
these applications introduce new and interesting properties of and extensions 
to BDEs. The climatic BDE model in Sect. 3, while keeping a small number of 
variables, introduces variables with more than two levels, as well as periodic 
forcing. Its solutions show that a simple BDE model can mimic rather well 
the solution set of a much more detailed model, based on nonlinear PDEs, 
as well as produce new and previously unsuspected results, such as a Devil's 
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Staircase and a "bizarre" attractor in phase-parameter space. 

The seismological BDE model in Sect. 4 introduces a much larger num- 
ber of variables, organized in a directed graph, as well as random forcing and 
state-dependent delays. This BDE model also reproduces a regime diagram 
of seismic sequences resembling observational data, as well as the results of 
much more detailed models [M] based on a large system of differential 
equations; furthermore it allows the exploration of seismic prediction meth- 
ods. 

3 A BDE Model for the El Nino/Southern 
Oscillation 

BDEs were first applied to paleoclimatic problems. Ghil et al. [1] used the 
exploratory power of BDEs to study the coupling of radiation balance of the 
Earth- atmosphere system, mass balance of continental ice sheets, and over- 
turning of the oceans' thermohaline circulation during glaciation cycles. On 
shorter time scales. Darby and Mysak [62] and Wohlleben and Weaver [63] 
studied the coupling of the sea ice with the atmosphere above and the ocean 
below in an interdecadal Arctic and North Atlantic climate cycle, respec- 
tively. Here we describe an application to tropical climate, on even shorter, 
seasonal-to-interannual time scales. 

The El-Nifio/Southern-Oscillation (ENSO) phenomenon is the most promi- 
nent signal of seasonal-to-interannual climate variability. It was known for 
centuries to fishermen along the west coast of South America, who witnessed 
a seemingly sporadic and abrupt warming of the cold, nutrient-rich waters 
that support the food chains in those regions; these warmings caused havoc 
to their fish harvests [Ml E5]. The common occurrence of such warming 
shortly after Christmas inspired them to name it El Nino, after the "Christ 
child." Starting in the 1970s, El Nino's climatic effects were found to be far 
broader than just its manifestations off the shores of Peru [Ml [66| [67] . This 
realization led to a global awareness of ENSO's significance, and an impetus 
to attempt and improve predictions of exceptionally strong El Nino events 
[681 [69]. 
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3.1 Conceptual ingredients 

The following conceptual elements are incorporated into the logical equations 
of our BDE model for ENSO variability. 

(i) The Bjerknes hypothesis: Bjerknes [70], who laid the founda- 
tion of modern ENSO research, suggested a positive feedback as a mechanism 
for the growth of an internal instability that could produce large positive 
anomalies of sea surface temperatures (SSTs) in the eastern Tropical Pa- 
cific. We use here the climatological meaning of the term anomaly, i.e., 
the difference between an instantaneous (or short-term average) value and 
the normal (or long-term mean). Using observations from the International 
Geophysical Year (1957-58), he realized that this mechanism must involve 
air-sea interaction in the tropics. The "chain reaction" starts with an initial 
warming of SSTs in the "cold tongue" that occupies the eastern part of the 
equatorial Pacific. This warming causes a weakening of the thermally direct 
Walker-cell circulation; this circulation involves air rising over the warmer 
SSTs near Indonesia and sinking over the colder SSTs near Peru. As the 
trade winds blowing from the east weaken and give way to westerly wind 
anomalies, the ensuing local changes in the ocean circulation encourage fur- 
ther SST increase. Thus the feedback loop is closed and further amplification 
of the instability is triggered. 

(ii) Delayed oceanic wave adjustments: Compensating for Bjerk- 
nes's positive feedback is a negative feedback in the system that allows a 
return to colder conditions in the basin's eastern part. During the peak of 
the cold-tongue warming, called the warm or El Nino phase of ENSO, west- 
erly wind anomalies prevail in the central part of the basin. As part of the 
ocean's adjustment to this atmospheric forcing, a Kelvin wave is set up in the 
tropical wave guide and carries a warming signal eastward. This signal deep- 
ens the eastern-basin thermocline, which separates the warmer, well-mixed 
surface waters from the colder waters below, and thus contributes to the 
positive feedback described above. Concurrently, slower Rossby waves prop- 
agate westward, and are reflected at the basin's western boundary, giving 
rise therewith to an eastward-propagating Kelvin wave that has a cooling, 
thermocline-shoaling effect. Over time, the arrival of this signal erodes the 
warm event, ultimately causing a switch to a cold. La Nina phase. 

( Hi ) Seasonal forcing: A growing body of work [5l [TH [72], [73], [71] [75] 
[76] points to resonances between the Pacific basin's intrinsic air-sea oscillator 
and the annual cycle as a possible cause for the tendency of warm events to 
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peak in boreal winter, as well as for ENSO's intriguing mix of temporal 
regularities and irregularities. The mechanisms by which this interaction 
takes place are numerous and intricate and their relative importance is not 
yet fully understood [76l [771 ESj • We assume therefore in the present BDE 
model that the climatological annual cycle provides for a seasonally varying 
potential of event amplification. 

3.2 Model variables and equations 

The model [6] operates with five Boolean variables. The discretization of 
continuous- valued SSTs and surface winds into four discrete levels is justified 
by the pronounced multimodality of associated signals (see Fig. lb of [H]). 

The state of the ocean is depicted by SST anomalies, expressed via a 
combination of two Boolean variables, Ti and T2. The relevant anomalous 
atmospheric conditions in the Equatorial Pacific basin are described by the 
variables Ui and 1/2- The latter express the state of the trade winds. For both 
the atmosphere and the ocean, the first variable, Ti or Ui, describes the sign 
of the anomaly, positive or negative, while the second one, T2 or U2, describes 
its amplitude, strong or weak. Thus, each one of the pairs (Ti, T2) and 
(f/i, U2) defines a four-level discrete variable that represents highly positive, 
slightly positive, slightly negative, and highly negative deviations from the 
climatological mean. The seasonal cycle 's external forcing is represented by 
a two-level Boolean variable S. 

The atmospheric variables Ui are "slaved" to the ocean [7H [79] : 

U,{t)=T,{t-P), 1 = 1,2. (15) 

The evolution of the sign Ti of the SST anomalies is modeled according to 
the following two sets of delayed interactions: 

(i) Extremely anomalous wind stress conditions are assumed to be nec- 
essary to generate a significant Rossby-wave signal R{t), which takes on the 
value 1 when wind conditions are extreme at the time and otherwise. By 
definition strong wind anomalies (either easterly or westerly) prevail when 
Ui = U2 and thus R{t) = Ui{t) A U2{t); here A is the binary Boolean opera- 
tor that takes on the value 1 if and only if both operands have the same value 
(see Sect. [2] and Table 1). A wave signal R{t) = 1 that is elicited at time t 
is assumed to re-enter the model system after a delay r, associated with the 
wave's travel time across the basin. Upon arrival of the signal in the eastern 
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equatorial Pacific at time t + r, the wave signal affects the thermocline-depth 
anomaly there and thus reverses the sign of SST anomalies represented by 



(ii) In the second set of circumstances, when R{t) = 0, and thus no signifi- 
cant wave signal is present, we assume that Ti(t+r) responds directly to local 
atmospheric conditions, after a delay /?, according to Bjerknes hypothesis; 
the delays associated with local coupled processes are taken all equal. 

The two mechanisms (i) and (ii) are combined to yield: 



here the symbols V and A represent the binary logical operators OR and 
AND, respectively (see Table 1). 

The seasonal-cycle forcing S is given by ^(t) = S(t — 1); the time t is 
thus measured in units of 1 year. The forcing 5* affects the SST anomalies' 
amplitude T2 through an enhancement of events when favorable seasonal 
conditions prevail: 



The model's principal parameters are the two delays (3 and r associated 
with local adjustment processes and with basin-wide processes, respectively. 
The changes in wind conditions are assumed to lag the SST variables by a 
short delay (3, of the order of days to weeks. For the length of the delay r we 
adopt Jin's [80] view of the delayed-oscillator mechanism and let it represent 
the time that elapses while combined processes of oceanic adjustment occur: 
it may vary from about one month in the fast-wave limit [HI [821 [83] to about 
two years. 

3.3 Model solutions 

Studying the ENSO phenomenon, we are primarily interested in the dynamics 
of the SST states, represented by the two- variable Boolean vector (Ti, T2). 
To be more specific, we deal with a four-level scalar variable 



Ti(t) = {[RAU^ ^t-r)}y {Rit-r)AU2{t-P)}; 



(16) 



T^it) = {[SATi] {t - /?)} V { [{SAT,) A T2] (t - /?)} 



(17) 



ENSO 



I 



—2, extreme La Nina, Ti = 0, T2 = 0, 

— 1, mild La Nifia, Ti = 0, T2 = 1, 

1, mild El Nifio, Ti = 1, T2 = 0, 

2, extreme El Nifio, Ti = 1, T2 = 1. 



(18) 



22 



In all our simulations, this variable takes on the values {—2, —1, 1, 2}, 
precisely in this order, thus simulating real ENSO cycles. The cycles follow 
the same sequence of states, although the residence time within each state 
changes as r changes at fixed (3. The period P of a simple oscillatory solution 
is defined as the time between the onset of two consecutive extreme warm 
events, ENSO = 2. We use the cycle period definition to classify different 
model solutions (see Figs. HHH]) • 

(i) Periodic solutions with a single cycle (simple period). Each 
succession of events, or internal cycle, is completely phase-locked here to the 
seasonal cycle, i.e., the warm events always peak at the same time of year. 
For each fixed (3, as r is increased, intervals where the solution has a simple 
period equal to 2, 3, 4, 5, 6 and 7 years arise consecutively. 

(ii) Periodic solutions with several cycles (complex period). 
We describe such sequences, in which several distinct cycles make up the full 
period, by the parameter P = P/n; here P is the length of the sequence and 
n is the number of cycles in the sequence. Notably, as we transition from a 
period of three years to a period of four years (see second inset of Fig. H]), P 
becomes a nondecreasing step function of r that takes only rational values, 
arranged on a Devil's Staircase. 

3.4 The quasi-periodic (QP) route to chaos in the BDE 
model 

The frequency- locking behavior observed for our BDE solutions above is a 
signature of the universal QP route to chaos. Its mathematical prototype is 
the Arnol'd circle map [H], given by the equation 

On+i = 9n + n + 2txK sin(27rft„) (mod 1). (19) 

Equation f[T^ describes the motion of a point, denoted by the angle 9 of its 
location on a unit circle, which undergoes fixed shifts by an angle Vt along 
the circle's circumference. The point is also subject to nonlinear sinusoidal 
"corrections," with the size of the nonlinearity controlled by a parameter K. 
The solutions of fll9p are characterized by their winding number 

u = u{Sl,K) = lim [(a„-^o)/n], 

n— >oo 

which can be described roughly as the average shift of the point per iteration. 
When the nonlinearity's influence is small, this average shift — and hence the 
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average period — is determined largely by fi; it may be rational or irrational, 
with the latter being more probable due to the irrationals' pervasiveness. As 
the nonlinearity K is increased, "Arnol'd tongues" — where the winding 
number uj locks to a constant rational over whole intervals — form and 
widen. At a critical parameter value, only rational winding numbers are left 
and a complete Devil's Staircase crystallizes. Beyond this value, chaos reigns 
as the system jumps irregularly between resonances [SH 

The average cycle length P defined for our ENSO system of BDEs is 
clearly analogous to the circle map's winding number, in both its definition 
and behavior. Note that the QP route to chaos depends in an essential way 
on two parameters: VL and K for the circle map and (3 and r in our BDE 
model. 

3.5 The "fractal sunburst": A "bizarre" attractor 

As the system undergoes the transition from an averaged period of two to 
three years a much more complex, and heretofore unsuspected, "fractal- 
sunburst" structure emerges (Fig. 0, and first inset in Fig. H]). As the wave 
delay r is increased, mini-ladders build up, collapse or descend only to start 
climbing up again. In the vicinity of a critical value (r = 0.5 years) that con- 
stitutes the pattern's focal point, these mini-ladders rapidly condense and 
the structure becomes self-similar, as each zoom reveals the pattern being 
repeated on a smaller scale. We call this a "bizarre" attractor because it is 
more than "strange": strange attractors occur in a system's phase space, for 
fixed parameter values, while this fractal sunburst appears in our model's 
phase-parameter space, like the Devil's Staircase. The structure in Fig. 4 
is attracting, though, only in phase space, for fixed parameter values; it is, 
therefore, a generalized attractor, and not just a bizarre one. 

The influence of the local-process delay /3, along with that of the wave- 
dynamics delay r, is shown in the three-dimensional "Devil's bleachers" (or 
"Devil's terrace," according to Jin et al. [71]) of Fig. [61 Note that the Jin et 
al. [751 [71] model is an intermediate model, in the terminology of modeling 
hierarchies [3], i.e. intermediate between the simplest "toy models" (BDEs 
or ODEs) and highly detailed models based on discretized systems of PDEs 
in three space dimensions, such as the general circulation models (GCMs) 
used in climate simulation. Specifically, the intermediate model of Jin and 
colleagues is based on a system of nonlinear PDEs in one space dimension 
(longitude along the equator). The Devil's bleachers in our BDE model 
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resemble fairly well those in the intermediate ENSO model of Jin et al. [71]. 
The latter, though, did not exhibit a fractal sunburst, which appears, on the 
whole, to be an entirely new addition to the catalog of fractals [86l Wi\ [88] . 

It would be interesting to find out whether such a bizarre attractor occurs 
in other types of dynamical systems. Its specific significance in the ENSO 
problem might be associated with the fact that a broad peak with a period 
between two and three years appears in many spectral analyses of SSTs and 
surface winds from the Tropical Pacific [89l[90]. Various characteristics of the 
Devil's Staircase have been well documented in both observations [90l[9T| [92] 
and GCM simulations [5l [73] of ENSO. It remains to see whether this will 
be the case for the fractal sunburst as well. 

4 A BDE Model for Seismicity 

Lattice models of systems of interacting elements are widely applied for mod- 
eling seismicity, starting from the pioneering works of Burridge and Knopoff 
P3] . Allegre et al. |H1], and Bak et al. [22]. The state of the art is summa- 
rized in [961 EH EH EHl [TOO]. Recently, colliding-cascade models [3 [HI [60l [6l] 
have been able to reproduce a wide set of observed characteristics of earth- 
quake dynamics jlOll I102[ I103j : (i) the seismic cycle; (ii) intermittency in 
the seismic regime; (iii) the size distribution of earthquakes, known as the 
Gutenberg-Richter relation; (iv) clustering of earthquakes in space and time; 
(v) long-range correlations in earthquake occurrence; and (vi) a variety of 
seismicity patterns premonitory to a strong earthquake. 

Introducing the BDE concept into the modeling of colliding cascades, 
we replace the elementary interactions of elements in the system by their 
integral effect, represented by the delayed switching between the distinct 
states of each element: unloaded or loaded, and intact or failed. In this 
way, we bypass the necessity of reconstructing the global behavior of the 
system from the numerous complex and diverse interactions that researchers 
are only mastering by and by and never completely. Zaliapin et al. [71 [8] 
have shown that this modeling framework does simplify the detailed study of 
the system's dynamics, while still capturing its essential features. Moreover, 
the BDE results provide additional insight into the system's range of possible 
behavior, as well as into its predictability. 
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4.1 Conceptual ingredients 

CoUiding-cascade models [H El EHl EI] synthesize three processes that play 
an important role in lithosphere dynamics, as well as in many other complex 
systems: (i) the system has a hierarchical structure; (ii) the system is con- 
tinuously loaded (or driven) by external sources; and (iii) the elements of the 
system fail (break down) under the load, causing redistribution of the load 
and strength throughout the system. Eventually the failed elements heal, 
thereby ensuring the continuous operation of the system. 

The load is applied at the top of the hierarchy and transferred downwards, 
thus forming a direct cascade of loading. Failures are initiated at the lowest 
level of the hierarchy, and gradually propagate upwards, thereby forming an 
inverse cascade of failures, which is followed by healing. The interaction of 
direct and inverse cascades establishes the dynamics of the system: loading 
triggers the failures, and failures redistribute and release the load. In its 
applications to seismicity, the model's hierarchical structure represents a fault 
network, loading imitates the effect of tectonic forces, and failures imitate 
earthquakes. 

4.2 Model structure and parameters 

(i) The model acts on a directed graph whose nodes, except the top one and 
the bottom ones, have connectivity six. Each node, except the bottom ones, 
is a parent to three children, that are siblings to each other. This graph is 
obtained from a directed ternary tree, which has its root in the top element, 
by connecting siblings, i.e., groups of three nodes that have the same parent. 

(ii) Each element possesses a certain degree of weakness or fatigue. An 
element fails when its weakness exceeds a certain threshold. 

(iii) The model runs in discrete time t = 0, 1, . . . . At each epoch a given 
element may be either intact or failed (broken), and either loaded or unloaded. 
The state of an element e at a epoch t is defined by two Boolean functions: 
Se{t) = 0, if an element is intact, and Se{t) = 1, if an element is failed; 
le{t) = 0, if an element is unloaded, and le(t) = 1, if an element is loaded. 

(iv) An element of the system may switch from one state {s, I) G {0, 1}^ 
to another under an impact from its nearest neighbors and external sources. 
The dynamics of the system is controlled by the time delays between the 
given impact and switching to another state. 

(v) At the start, t = 0, all elements are in the state (0, 0), intact and un- 
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loaded. Most of the changes in the state of an element occur in the following 
cycle: 

(0,0)^(0,1)^(1,1)^(1,0)^(0,0)... 

Other sequences, however, are also possible, except that a failed and loaded 
element may switch only to a failed and unloaded state, (1, 1) — »■ (1,0). The 
latter transition mimics fast stress drop after a failure. 

(vi) All the interactions take finite, nonzero time. We model this by 
introducing four basic time delays: A^,, between being impacted by the load 
and switching to the loaded state, (-,0) {'A)] ^f, between the increase 
in weakness and switching to the failed state, (0,-) — (1,-); Ao, between 
failure and switching to the unloaded state, (-, 1) (-,0); and Ah, between 
the moment when healing conditions are established and switching to the 
intact (healed) state, (1, ■) — ^ (0, ■). 

The duration of each particular delay, from one switch of an element's 
state to the next, is determined from these basic delays, depending on the 
state of the element as well as of its nearest neighbors during the preceding 
time interval (see [7] for details). This represents yet another generalization 
of the set of deterministic, autonomous equations with fixed delays Oij: 
here the effective delays are both variable and state-dependent. 

(vii) Failures are initiated randomly within the elements at the lowest 
level. 

The two primary delays in this system are the loading time A^, necessary 
for an unloaded element to become loaded under the impact of its parent, 
and the healing time necessary for a broken element to recover. 

Conservation law. The model is forced and dissipative, if we associate 
the loading with an energy influx. The energy dissipates only at the lowest 
level, where it is transferred downwards, out of the model. In any part of 
the model not including the lowest level, energy conservation holds, but only 
after averaging over sufficiently large time intervals. On small intervals it 
may not hold, due to the discrete time delays involved in energy transfer. 

Model solutions. The output of the model is a catalog C of earthquakes 
— i.e., of failures of its elements — similar to the simplest routine catalogs 
of observed earthquakes: 

C = {tk, nik, hk), k = 1,2,...; 4 < tk+i. (20) 

In real-life catalogs, tk is the starting time of the rupture; rrik is the magni- 
tude, a logarithmic measure of energy released by the earthquake; and hk is 
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the vector that comprises the coordinates of the hypocenter. The latter is 
a point approximation of the area where the rupture started. In our BDE 
model, earthquakes correspond to failed elements, rrik is the level at which 
the failed element is situated within the directed graph, while the position 
of this element within its level is a counterpart of hk- 

4.3 Seismic regimes 

A long-term pattern of seismicity within a given region is usually called a seis- 
mic regime. It is characterized by the frequency and irregularity of the strong 
earthquakes' occurrence, more specifically by (i) the Gutenberg-Richter re- 
lation, i.e. the time-and-space averaged magnitude-frequency distribution; 
(ii) the variability of this relation with time; and (iii) the maximal possi- 
ble magnitude. The notion of seismic regime here is a much more complete 
description of seismic activity than the "level of seismicity," often used to dis- 
criminate among regions with high, medium, low and negligible seismicity; 
the latter are called aseismic regions. 

The seismic regime is to a large extent determined by the neotectonics of 
a region; this involves, roughly speaking, two factors: (i) the rate of crustal 
deformations; and (ii) the crustal consolidation, determining what part of 
deformations is realized through the earthquakes. However, as is typical for 
complex processes, the long-term patterns of seismicity may switch from one 
to another in the same region, as well as migrate from one area to another 
on a regional or global scale [104[ 1105] . Our BDE model produces syn- 
thetic sequences that can be divided into three seismic regimes, illustrated 
in Figs. EHm 

Regime H: High and nearly periodic seismicity (top panel of Figs. [7] and 
IH]). The fractures within each cycle reach the top level, m = L, where our 
underlying ternary graph has depth L = 6. The sequence is approximately 
periodic, in the statistical sense of cyclo-stationarity |106j . 

Regime I: Intermittent seismicity (middle panel of Figs. [7] and [8]). The 
seismicity reaches the top level for some but not all cycles, and cycle length 
is very irregular. 

Regime L: Medium or low seismicity (lower panel of Figs. [7] and [8]). No 
cycle reaches the top level and seismic activity is much more constant at a 
low or medium level, without the long quiescent intervals present in Regimes 
H and I. 

The location of these three regimes in the plane of the two key parameters 
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(A^, Ah) is shown in Fig. [101 Figures [7HT^ were computed for a tree depth 
of L = 6, i.e. 1093 nodes. Many calculations were also carried out for L = 7, 
i.e. 3280 nodes, and the results were similar, but are not reported here. 



4.4 Quantitative analysis of regimes 

The quantitative analysis of model earthquake sequences and regimes is fa- 
cilitated by the two measures described below. 

Density of failed elements. The density p{t) of the elements that are 
in a failed state at the epoch n is given by: 



Here I'mit) is the fraction of failed elements at the m-th level of the hierarchy 
at the epoch t, while L is the depth of the underlying tree. Sometimes we 
consider this measure averaged over a time interval, or a union of intervals, 
/ and denote it by p{I) . The density p{t) for the three sequences of Fig. [7] is 
shown in Fig. [HI 

Irregularity of energy release. The second measure is the irregularity 
of energy release over the time interval /. It is motivated by the fact 
that one of the major differences between regimes resides in the temporal 
character of seismic energy release. The measure Q is defined by the following 
sequence of steps: 

(i) First, define a measure of seismic activity within the time interval, 
or union of time intervals, / as 



The summation in (122!) is taken over all events within /, i.e., ti G /; nj is the 
total number of such events, and rrii is the magnitude of the i-th event. The 
value of B equalizes, on average, the contribution of earthquakes with dif- 
ferent magnitudes, that is from different levels of the hierarchy. In observed 
seismicity, has a transparent physical meaning: given an appropriate 

choice of B, it estimates the total area of the faults unlocked by the earth- 
quakes during the interval / |107j . This measure is successfully used in several 
earthquake prediction algorithms [96] . 

(ii) Consider a subdivision of the interval / into a set of nonoverlapping 
intervals of equal length e > 0. For simplicity we choose e such that |/| = eiV/, 



p{t) = [iy,{t) + --- + ULit)]/L. 



(21) 




i=l 



(22) 
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where | ■ | denotes the length of an interval and Nj is an integer. Therefore, 
we have the following representation: 

Ni 

I=[jlj, |4| = e, k = l,...,Nj; /,n4 = 0forjVfc- (23) 
(iii) For each k = 1, . . . , Nj we choose a fc-subset 

i=il,...,ik 

that maximizes the value of the accumulated S: 



j:[Q{k)] = J:*{k) = max <^ S 



(24) 



Here the maximum is taken over all /c-subsets of the covering set fl23l) . 
(iv) Introducing the notations 

f:{k) = S*(A;)/S(J), T{k) = ke/\I\, 



(25) 



we finally define the measure Q of clustering within the interval / as 

^(/) = ^max^^{S(fc)-r(A;)}. (26) 

Figure M illustrates this definition by displaying the curves S — r vs. r 
for the three synthetic sequences shown in Fig. [71 The curves give, essen- 
tially, the maximum seismic activity minus the mean activity, as a function 
of length of time over which the activity occurs, and the maximum of each 
curve gives the corresponding value of Q. The more clustered the sequence, 
the more convex is the corresponding curve, the larger the corresponding 
value of Q, and the shorter the interval for which this value of Q is realized. 
Despite its somewhat elaborate definition, Q has a transparent intuitive in- 
terpretation: it equals unity for a catalog consisting of a single event (delta 
function, burst of energy), and it is zero for a marked Poisson process (uni- 
form energy release). Generally, it takes values between and 1 depending 
on the irregularity of the observed energy release. 
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4.5 Bifurcation diagram 

Figure [TT] provides a closer look at the regime diagram of Fig. (TO) it illustrates 
the transition between regimes in the parameter plane {Al,Ah)- To do so, 
Fig. [11] (a) shows a rectangular path in the parameter plane that passes 
through all three regimes and touches the triple point. We single out 30 
points along this path; they are indicated by small circles in the figure. The 
three pairs of points that correspond to the transitions between regimes are 
distinguished by larger circles and marked in addition by letters, for example 
(A) and (B) mark the transition from Regime H to Regime L. 

We estimate the clustering G{I) and average density p{I) over the time in- 
terval / of length 2-10^ time units, for representative synthetic sequences that 
correspond to the 30 marked points along the rectangular path in Fig. [TTh . 
Figure [TTb is a plot of p(/) vs. for these 30 sequences. The values 

of Q drop dramatically, from 0.8 to 0.18, between points (A) and (B): this 
means that the energy release switches from highly irregular to almost uni- 
form between Regimes H and L. This transition, however, barely changes 
the average density p of failures. 

The transitions between the other pairs of regimes are much smoother. 
The clustering drops further, from Q = 0.18 to ^ ^ 0.1, and then remains 
at the latter low level within Regime L. It increases gradually, albeit not 
monotonically, from 0.1 to 0.8 between points (C) and (A), on its way through 
regimes I and H. The increase of along the right side of the rectangular 
path in Fig. [TTb . between points (F) and (A), corresponds to a decrease of p 
and a slight increase of clustering Q, from 0.5-0.6 to ^ 0.8. 

The transition between regimes is illustrated further in Fig. [121 Each 
panel shows a fragment of the six synthetic sequences that correspond to 
the points (A)-(F) in Fig. [TTk . The sharp difference in the character of 
the energy release at the transition between Regimes H (point (A)) and L 
(point (B)) is very clear, here too. The other two transitions, from (C) to 
(D) and (E) to (F), are much smoother. Still, they highlight the intermittent 
character of Regime I, to which points (D) and (E) belong. 

Zaliapin et al. [S] considered applications of these results to earthquake 
prediction. These authors used the simulated catalogs to study in greater 
detail the performance of pattern recognition methods tested already on ob- 
served catalogs and other models [961 [M [M [Ml EOl nH [ml [ITSl El] , 
devised new methods, and experimented with combination of different indi- 
vidual premonitory patterns into a collective prediction algorithm. 
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5 BDEs on a Lattice and Cellular Automata 
(CAs) 



While the development and applications of BDEs started about 25 ago, this is 
a very short time span compared to the long history of ODEs, PDEs, maps, 
and even CAs. The BDE results obtained so far, though, are sufficiently 
intriguing to warrant further exploration. In this section, we provide some 
preliminary results on BDE systems with a large or infinite number of 
variables, and we discuss in greater detail their connections with CAs fU] 
[IHlIll, (see also Fig. 1 and Sect. 2.6). 

These "partial BDEs" that we are led to explore, with N oo Boolean 
variables, were mentioned in passing in [2j, and stand in the same relation to 
"ordinary BDEs," explored so far, as PDEs do to ODEs. The classification 
of what we could call now ordinary BDEs into conservative and dissipative 
(Sect. [2]) suggests that partial BDEs of different types can exist as well. 



5.1 Towards partial BDEs of hyperbolic and parabolic 
type 

We want, first of all, to clarify what the "correct" BDE equivalent of partial 
derivatives may be. To do so, we start by studying possible candidates for 
hyperbolic and parabolic partial BDEs. Intuitively, these should correspond 
to generalizations of conservative and dissipative BDEs, respectively (see 
Sects. 2.3 and 2.5), that are infinite-dimensional (in phase space). We are 
thus looking for the discrete-variable version of the typical hyperbolic and 
parabolic PDEs: 

g-^^i^^t) = ^^(^'^)' (27) 

^^(^'^) = ^^(^'^)' (28) 

where, in the spatially one-dimensional case, v : M."^ ^ M. We wish to 
replace the real- valued function v with a Boolean function m : — B, i.e. 
u{z,t) = 0, 1, with 2; G M and t G M"^, while conserving the same qualitative 
behavior of the solutions. 

Since for Boolean-valued variables \x—y\ = x\jy (see Sect.[2]and Table 1), 
one is tempted to use the "eXclusive OR" (XOR) operator y for evaluating 
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differences. Moreover, one is led to introduce a time delay 6t and a space 
delay 6^ when approximating the derivatives in Eqs. (1271) and (125]) by finite 
differences. First-order expansions then lead to the equations: 

u{z, t + et)\/ u{z, t) = u{z + e,, t) V u{z, t) (29) 

in the hyperbolic case, and 

u{z, t + et)v u{z, t) = u{z - e,, t) V u{z + e,, t) (30) 

in the parabolic one. 

5.2 Boundary conditions and discretizing space 

The pure Cauchy problem for Eq. fl2^ on the entire real line [115j . z G 
(—00, 00) has solutions of the form: 

u{z,t) =u{z + 9„t-et). (31) 

This first step towards a partial BDE equivalent of a hyperbolic PDE displays 
therewith the expected behavior of a "wave" propagating in the {z, t) plane. 
The propagation is from right to left for increasing times when the "plus" sign 
is chosen in the right-hand side (rhs) of Eq. (1291) . as we did, but it could be in 
the opposite direction choosing the "minus" sign instead. The solution (1311) 
of fl29l) exists for all times t > 9t and it is unique for all delays {9^, 9t) G (0, 1]^, 
and for all initial data Uq{z, t) with z G (— cxd, oo) and t G [0, 9t). 

In continuous space and time, Eq. fl29|) under consideration is conserva- 
tive and invertible. Aside from the pure Cauchy problem discussed above, 
when UQ^z.t) is given for z G (—00,00), one can also formulate for fl2I?l) a 
space-periodic initial boundary value problem (IBVP), with uq{z^ t) given for 
z G [0,Tj) and uq{z + Tz,t) = uo{z,t). The solution of this IBVP displays 
periodicity in time as well, with Tt = Tz9t/9z. This time-and-space periodic 
solution exists for all time and is unique under conditions that are analogous 
to those stated for the pure Cauchy problem above. 

Next, we analyze a discrete version of Eq. (I29p . which is obtained by 
studying the evolution of the system on a 1-D lattice. Specifically, one con- 
siders the grid {zi = zo + i9z}i(zz for a fixed 2:0 G M and assumes that the 
initial state is constant within the space intervals Jj = [zq + i9z, ZQ + {i + 1)9^] 
for t G [0,^(), i.e., Uo{z,t) = %,i(i)Ii, where Ij is the characteristic 
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function of the interval equal to unity on Jj and to zero outside it. Cor- 
respondingly, the behavior of u{z,t) is determined by the evolution of the 
elements of the set {ui(t) = u{zi,t)}i(zz and one gets: 

Ui{t) = Ui^i{t-et), ieZ. (32) 

Equation fl32l) is a simple linear BDE system ([7]), with Cij = and equal 
delays 6ij = 6t for all i and j, in the limit n oo. The discretization thus 
makes it more evident that, in the IBVP case, where {no,i(t) = noi+Ar(t)}igz, 
the solution is immediately periodic in time, without transients, and with 
period Tt = NOt- Notice that for all {Oz,0t) G (0,1]^, whether rational or 
not, one can choose the space discretization generated by multiples of 6z and 
the resulting partial BDE system with constant initial data within the space 
intervals will depend on the single delay 6t. The same observation should 
apply to more general cases, and therefore partial BDEs, whether hyperbolic 
or not, but containing only one space delay 6z and one time delay 6t cannot 
display solutions with increasing complexity. The case of partial BDEs with 
higher time derivatives and admitting therewith more than one time delay, 
is left for future work. 

We merely note here that approximating in the hyperbolic PDE fl27I) 
to the second order |116] yields the same partial BDE fl5Ul) that was obtained 
from the first-order approximation of the parabolic PDE: 

uiz, t) = u{z - Bz. t-Qt)\/ u{z, t-et)y u{z + ez,t- Ot). (33) 

The equivalence is apparent by using the associativity of addition (mod 2) 
in B (see Sect. [2]). 

This result is slightly, but not utterly unexpected: it is fairly well known 
in the numerical analysis of PDEs (i) that higher-order finite-difference ap- 
proximations of partial derivatives can lead to PAEs (see Fig. 1) that are 
consistent with a different PDE, containing higher-order derivatives [117] : 
and (ii) that such higher-order approximations, when they are stable, are 
more likely than not to be dissipative. Still, Eqs. fl29|) . fl30l) and fl33|) show 
that finding the "correct" partial BDE equivalent of a given PDE is not quite 
trivial. To get further insight into the behavior of Eq. (1551) . or, equivalently, 
Eq. flHU]) . we consider again the case of discrete space, and show that this 
system is equivalent, in turn, to a particular CA in the limit of infinite size. 
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5.3 Partial BDE (BSD and EGA rule 150 



Applying the space discretization scheme from the previous subsection (Sect. [5^ 
and using the same space-periodic initial data, one finds that 

Ui{t) =Ui_i{t-9t)\/u,{t-9t)s^u,+i{t-9t), leZ. (34) 

In order to establish the equivalence of Eq. flMl) with the general form f[T^ of 
EGA evolution (in one space dimension), one needs to specify the boundary 
conditions. For periodic boundary conditions Ui = Ui+N, and the number n of 
BDEs that we can associate with such an EGA is n = 2N + 1; this number, 
though, can be slightly lower or higher if one uses Dirichlet or Neumann 
boundary conditions. 

To identify the particular Boolean function that gives the EGA corre- 
sponding to Eq. flM|l . recall that in an EGA, i.e. a 1-D GA where the 
interactions involve only the nearest right and left neighbors Ui±i, the sin- 
gle rule valid at all sites can be described by a binary string. This string 
summarizes the truth table of the rule, by assigning the values of the inputs 
(wj-i, Mj, Mj+i) in decreasing order, from 111 to 000. Gorrespondingly, one 
gets 2^ = 8 binary digits, 1 or 0, for each possible output. The 8-digit string 
that characterizes the 3-site rule on the rhs of Eq. (!34l) is 10010110. For 
brevity. Wolfram replaces [11] each such 8-digit binary number by its deci- 
mal representation, which yields, in this case, rule 150 [ID]. We also recall 
that, for a finite number of variables, i.e. for i G [— A^, N] with finite size A^, 
the EGA version of our pigeon-hole lemma (Theorem 2.3 in Sect. 2.2) states 
that all solutions of such an automaton have to become stationary or purely 
periodic after a finite number of time steps. For A^ infinite, however, rule 
150 yields interesting behavior, with self-similar, fractal patterns embedded 
in its spatio-temporal structure. 

Next, we study in detail the simple case of Eq. (134|) with the time-constant 
initial state, Mo,i(t) = Mi(0) for < t < 9t, and assume 9t = 1 without loss 
of generality, and 9z = 9t for simplicity. We thus verify merely that our 
partial BDE does generate the rule-150 EGA behavior, which is already quite 
interesting, while we expect less trivial BDEs to yield distinct, and possibly 
even more interesting, behavior. 

We notice that the system's dynamics for any initial state of the pure 
Gauchy problem can be obtained from the evolution of the corresponding 
EGA that starts from initial data with a single nonzero value, which we 
show in Fig. [T5h . This property of solutions is due to the fact that rule 150 
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exhibits the important simphfying feature of "additive superposition" [ID], 
as evident from Fig. [T5b . where the colhsion of two "waves" is plotted; this 
feature is the resuh of hnearity (mod 2), as discussed in [21 [3] and in Sect. [2] 
here. 

In the IBVP case, the solutions can behave quite differently, depending 
on the value of and on the initial state. We report here only the results for 
the IBVP with periodic initial data Ui{0) = Ui^xA^) number of variables 
2N + 1, with chosen to be an integer multiple of the space period T^. In 
this case — apart from the presence of possible fixed points, in particular for 
Tz < 3 — the solutions of Eq. (IMI) display the longest time periods when 
Tz = N. More precisely, Mo(0) = mt-(O) = 1 and Ui{0) = for i not a multiple 
of Tz] this choice of the initial state has the longest-possible space period at 
fixed size N. 

In agreement with other known results for EGA 150 [iQl I118j . we find 
that our solutions with = N are immediately time periodic for A^ not 
a multiple of three, whereas there is a transient when A^ = 3p,p G N; this 
transient is of length 1 for A^ odd and of length 2^~^, where 2^ is the largest 
power of two which divides A^, otherwise. The EGA with rule 150 belongs 
to the third class in Wolfram's classification [iQl HI], in which evolution for 
N ^ oo can lead to aperiodic, chaotic patterns. While the behavior of this 
EGA is predictable from the knowledge of the initial state for the IBVP with 
finite A^, the length of the time period can increase rapidly with A^. In 
particular, one finds a time period already as long as 511 for space-periodic 
initial data with = N = 19. 

To show that the behavior of the IBVP for Eq. fl33l) is not periodic in time 
in the limit of A^ ^ oo, we present in Fig. [TUthe results for a random initial 
state of length A^ = 100. These results do not show any "recurrent pattern," 
apart from the expected [TTl |13] appearance of characteristic "triangles" that 
still emerge from the chaotic distribution of empty and occupied sites. 

We have seen in Sect. 2 (see Definition 2.1) that ordinary BDEs are 
conservative if their behavior is immediately periodic for any initial data, 
and that correspondingly they are also time-reversible. Our findings show 
that Eq. (155]) is dissipative, since one observes transients for the periodic 
IBVP with various A^ values; this is expected from its equivalence with the 
"parabolic" BDE (|30|1 . Nevertheless, in the ordinary BDE context, the con- 
nectives of the form (l33l) . based on the y operator, are often conservative at 
least for some set of distinct delays. We wish therefore to extend the analysis 
of this partial BDE — from the case of all delays being equal, 9ij = Of, as 
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above — to the case of distinct 9t values in the rhs of Eq. 0341) . since this could 
lead to conservative behavior and to increasingly complex solutions for some 
set of delays. Doing so, however, goes beyond the purpose of the present 
paper and will have to await future work. 

Another attempt at formulating parabolic partial BDEs is given by the 
equations: 

u{z,t) = [u{z-e,,t-et)yu{z + e,,t-et)]yu{z,t-et), (35) 

u{z,t) = [u{z-e,,t-9t) Au{z + e,,t-9t)]\7uiz,t-et), (36) 

obtained from the parabolic PDE (l28l) by replacing dzz with the OR and the 
AND operator, respectively. From the known results on ordinary BDEs, the 
solutions in these cases can be expected to reproduce closely the behavior 
of dissipative PDEs. This intuitive conjecture is confirmed by the analysis 
in terms of their EGA equivalent. Thus, Eq. fl35|) corresponds to rule 54, in 
the first class, and Eq. fl36|) to rule 108, in the second class. Accordingly, 
the evolution of their solutions leads to fixed points or to small-period limit 
cycles, respectively; in neither case does the N —>■ oo limit display chaotic 
behavior. We thus expect these connectives to be good candidates as simple 
examples of "correct" partial BDE equivalents of dissipative PDEs. 

5.4 Summary on partial BDEs and future work 

Our results on classifying finite BDE systems, on the one hand, and our 
replacing partial derivatives in PDEs by Boolean operators, on the other, 
seem to provide interesting insights into the correspondence between partial 
BDEs and PDEs. In Sects. 5.1-5.3 we have only considered relatively easy 
cases, in which close correspondence exists between our partial BDEs and 
EGAs. This correspondence sheds new light on the known results of certain 
cellular automata. 

We summarize our results on partial BDEs in Table 2. In the degenerate 
case of all delays being equal, 6ij = 6t for all i and j, partial BDEs are always 
equivalent to particular GAs in the limit of infinite size, and their expected 
behavior can be obtained from the analogy. We gave examples of dissipative 
partial BDEs corresponding to EGAs in different classes of Wolfram's classifi- 
cation [11] . In particular, a first-order approximation of the spatial derivative 
in the parabolic PDE fl28|) . as well as a second-order approximation of this 
derivative in the hyperbolic PDE fl27|) yielded the same partial BDE fl33|) . 
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equivalent to EGA 150, which can behave very differently depending on the 
starting data. On the other hand, Eq. fl29p gave immediately periodic so- 
lutions for all the starting data; hence this partial BDE is conservative. In 
the binary string formulation reviewed in Sect. 15. 3[ it corresponds to EGA 
rule 170, which is among the "forbidden" rules of Wolfram [11] because it 
describes asymmetric interactions. Given its simplicity, though, Eq. fl2^ and 
hence rule 170 seem perfectly legitimate in the present context. 

More difficult situations, such as the case of nonconstant starting data in 
the initial interval, are left for future studies. To give an idea of the possible 
outcomes, consider the solution of Eq. (13T|) with initial data of the Riemann 
type, with a single jump in Uo{z, t) at z = and to = 1/2 < 9t, i.e., Uq{z, t) = 
1 for — oo < 2; < and < t < 1/2, and UQ{z,t) = in the rest of the ini- 
tial strip, {{z,t) : < z < 00, t < et}[j {{z, t) : -00 < z <0, 1/2 <t < 9t}. 
An obvious conjecture is that the solution will still be a right- or left-traveling 
wave, depending on the sign taken in the rhs of the partial BDE. More gen- 
erally, certain ^-periodic IBVPs will also be well posed and, in the more 
intriguing case of Eq. (1551) . additive superposition will provide insight into 
solution behavior, which might include certain solutions that exhibit both z- 
and t-periodicity. 

In the parabolic partial BDE (l35l) we conjecture instead that the solution 
u = 1 will be asymptotically stable, even for rationally independent 6^ and 
6t and nonconstant initial data, at least for large \z\. Gonversely, u = 
should be asymptotically stable when replacing the OR operator by AND 
as in Eq. fl36|l : see Theorem 2.9 in Sect. 2.5. However, nontrivial solutions 
of parabolic problems could be obtained in the presence of time-constant or 
time-periodic forcing, like for PDEs. 

One further step would consist of looking at the same connectives with 
different delays, i.e. different Of values in the variables in the rhs of the equa- 
tions. This should in particular allow one to better understand the classifica- 
tion of partial BDEs into conservative and dissipative, possibly depending on 
the open set of delays under consideration. We conjecture, moreover, that 
— ^for distinct, irrationally related delays — connectives similar to the one 
in Eq. fl33l) could have solutions of increasing complexity. If so, this would 
provide us with an even richer metaphor for evolution than either ordinary 
BDEs or EGAs. 
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6 What Next? 



The most promising development in the theory and apphcations of BDEs 
seems to be the extension to an infinite, or very large, number of variables 
as discussed in Sect. 4 [n ^ 10^) and Sect. 5 {n ^ oo). Inhomogeneous 
partial BDEs can be easily handled by introducing "variable coefficients," 
i.e. multiplication of the right-hand side by site-dependent Boolean variables 
or functions. 

Another important possibility is the randomization of various aspects of 
the BDE formulation. Wright et al. [119] have already considered ensemble 
averaging over BDE solutions with randomized initial data, while Zaliapin et 
al. [3 [S] have considered random forcing. It would be even more interesting 
to consider the random perturbation of delays, first in ordinary and then in 
partial BDEs; see also work along these lines in Boolean networks |58j . 

Next, an important but harder goal of the theory will be to develop in- 
verse methods. In other words, given the behavior of a complicated natural 
system, which can be described in a first approximation by a finite number 
of discrete variables, one would like to discover a good connective linking 
these variables and yielding the observed behavior for plausible values of the 
delays. In this generality, the inverse problem for ODEs, say, is clearly in- 
tractable. But some of the direct results obtained so far — see, for instance, 
the asymptotic-simplification results in Sect. 2.5 — hold promise for BDEs, 
at least within certain classes and with certain additional conditions. In- 
tuitively, the behavior of BDEs, although surprisingly rich, is more rigidly 
constrained than that of fiows. Certain inverse-modeling successes have also 
been reported for cellular automata; see [TT] and references there. 

From the point of view of applications, BDEs have been applied fairly 
extensively by now to climate dynamics [H [62l [631 1119i I120[ I121[ I122j and are 
making significant inroads into solid-earth geophysics [Tj, [8] • Most interesting 
are recent applications to the life sciences (Neumann and Weisbuch |123[I124] . 
Gagneur and Casari [13], Oktem et al. [59]), which represent in a sense a 
return to the motivation of the geneticist Rene Thomas, originator of kinetic 
logic [201 [Ml EH]; see also Sect. 2.6 and Appendix A for details. 

BDEs may be suited for the exploration of poorly understood phenomena 
in the socio-economic realm as well. Moreover, the robustness of fairly regular 
solutions in a wide class of BDEs, for many sets of delays and a variety of 
initial states, suggests interesting applications to certain issues in massively 
parallel computations. 
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Appendix A: BDEs and Kinetic Logic 

A mathematical model closely related to BDEs was formulated by R. Thomas 
[3S] in a slightly different form. The difference can best be explained in 
terms of "memorization variables" Xij{t) = Xj{t — 6ij). In our formulation 
Xij{t) is a purely delayed state variable; that is, Xij{t) is determined only by 
the state of the system at time t — 6ij. Thomas allows the memorization 
variable Xij{t) to depend on the state of the system up to and including time 
t: when a change takes place in any variable x^ in the time interval {t — 6ij, t), 
then the memorization variable Xij may change as well. 

Specifically, if for some t' E {t — 9ij, t) the state of the system is such that 

Xjit') = fjMt'), . . .,Xnit')) ^ x,{t - %), (37) 

then the memorization variable Xij{t), previously equal to Xj(t — %), is 
changed to 

x.jit) =Xj{t'). (38) 

The adjustment of the variables Xij is in effect selectively erasing some of 
the memory of the system. The resulting solutions are usually simpler and 
hence easier to study than the typical solution of our BDEs. Referring to 
the examples shown in Figs. 2 and 3 here, the increasing complexity of the 
solution reflects the fact that the memory of the system contains more and 
more information as time goes on. 
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In Thomas's kinetic-logic formulation, such solutions of increasing com- 
plexity cannot arise. They are due precisely to the "conflicts" between vari- 
able values that were avoided on purpose in kinetic logic [13]. Some recent 
work on Boolean networks [581159] use a similar policy of eliminating behavior 
that appears to be too complicated, at least for certain purposes, and thus 
too hard to capture numerically. Eliminating the selective memory erasure 
of Eqs. (jnZD and (1551) seems, on the other hand, to provide a cleaner, richer 
and more versatile mathematical theory [Hill [3]. 
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Table 1. Common Boolean operators 



Mathematical 


Engineering 


Name 


Description 


symbol 


symbol 






X 


NOT X 


negator 


not true when x is true 


x\/ y 


X OR y 


logical or 


true when either x or y ot both are true 


X Ay 


X AND y 


logical and 


X Ay = (x V y) 


x\/y 


X XOR y 


exclusive or 


true only when x and y are not equal 


X Ay 






true only when x and y are equal 
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Table 2. Results on partial BDEs 



PDE 


Approxim- 


PBDE 


EGA 


EGA 


Eq. and 


dtv = 


ation 


Uiit + Ot) = 


rule 


class 


behavior 


dzV 


V (1st order) 




170 




(El G 


dzzV 


V 




54 


1 


(ESD D 


dzzV 


A 


[n,_i(t) A Ui+i{t)] V Ui{t) 


108 


II 


iSD D 


dzzV 

dzV 


V (1st order) 

V (2nd order) 


[ui_i{t) V Ui+i{t)] V Ui{t) 


150 


III 


(I33D D 



Summary of results on the partial BDEs (PBDEs) obtained from differ- 
ent approximations for the spatial derivative in the simplest parabolic and 
hyperbolic PDEs. The temporal derivative is always approximated to first 
order by the V operator. The last column gives the equation number and the 
behavior of the solutions, conservative (G) or dissipative (D). Notice that, 
though all but Eq. (!29|) in this table are dissipative, it is only Eq. (133|) that 
displays chaotic behavior in the limit of infinite lattice size. 
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Automata 

X discrete, t discrete 
, {Turing machines, real computers, 
CAs, conservative logic) 



Figure 1: Tlie place of BDEs within dynamical system theory. Note the 
links: The discretization of t can be achieved by the Poincare map (P-map) 
or a time-one map, leading from Flows to Maps. The opposite connection 
is achieved by suspension. To go from Maps to Automata we use the 
discretization of x. Interpolation and smoothing can lead in the opposite 
direction. Similar connections lead from BDEs to Automata and to Flows, 
respectively. Modified after Mullhaupt [2]. 
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Figure 2: Solutions of the system of two conservative BDEs ([6]) for the delay 
6 = 0.977 and < t < 40 [Tj. The tick marks on the t-axis indicate the 
times at which jumps in either xi or X2 take place. After Dee and Ghil [1]. 
Copyright ©Society for Industrial and Applied Mathematics; reprinted with 
permission. 
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Figure 3: Jump function J{t) for the particular solution of a conservative 
BDE with the irrational delay 9 = (-\/5 — 1)/2; see Eq. ([9]) and further details 
in the text. Reproduced from [3], with kind permission of Springer Science 
and Business Media. 



57 




—I i 1—" 

1.15 1.20 1.25 
Wave delay T (yrs) 



1 I I I I I I I I I I I 

0.2 0.6 1.0 1.4 1.8 2.2 
Wave delay i (yrs) 

Figure 4: Devil's staircase and fractal sunburst for a BDE model of the El- 
Nino/Southern-Oscillation (ENSO) phenomenon. Plotted in the bifurcation 
diagram is the average cycle length P vs. the wave delay r for a fixed local 
delay /3 = 0.17. Blue dots indicate purely periodic solutions; orange dots are 
for complex periodic solutions; small black dots denote aperiodic solutions. 
The two insets show a blow-up of the overall, approximate behavior between 
periodicities of two and three years ( "fractal sunburst" ) and of three and four 
years ("Devil's staircase"). Modified after Saunders and Ghil [6]. 
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Figure 5: Fractal sunburst: a BDE solution pattern in phase-parameter 
space, for a dissipative BDE system with periodic forcing. The plot is a blow- 
up of the transition zone from average periodicity two to three years in Fig. HI 
here r = 0.44-0.58,/? = 0.17. The inset is a zoom on 0.490 < r < 0.504. 
A complex mini- staircase structure reveals self-similar features, with a focal 
point at r ^ 0.5. Modified after Saunders and Ghil [6]. 
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Figure 6: The Devil's bleachers in our BDE model of ENSO. The three- 
dimensional regime diagram shows the average cycle length P, portrayed in 
both height and color, vs. the two delays (3 and r. Oscillations are produced 
even for very small values of (3, as long as /3 < r. Variations in r determine 
the oscillation's period, while changing (5 establishes the bottom step of the 
staircase, shifts the location of the steps, and determines their width. After 
Saunders and Ghil [6]. 
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Figure 7: Three seismic regimes in a coUiding-cascade model of lithospheric 
dynamics; each earthquake sequence illustrates characteristic features of the 
corresponding regime and only a small fraction of each sequence is shown. 
Top panel - regime H (High), Ah = 0.5 ■ 10^; middle panel - regime I 
(Intermittent), Ah = 10^; bottom panel - regime L (Low), Ah = 0.5 ■ 10^; 
in all three panels Al = 0.5 ■ 10^ (see also Fig. [TOl below) and the number of 
nodes in the simulated lattice is 1093, for a tree depth of L = 6, the maximum 
magnitude of any earthquake in the BDE model. Reproduced from Zaliapin 
et al. [7j, with kind permission of Springer Science and Business Media. 
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Figure 8: Three seismic regimes: internal dynamics of the BDE model. The 
panels show the density p{n) of broken elements in the system, as defined by 
Eq. (12T!) : they correspond to the synthetic sequences shown in Fig. [3 Top 
panel - Regime H; middle panel - Regime I; and bottom panel - regime L. 
Reproduced from Zaliapin et al. [7] , with kind permission of Springer Science 
and Business Media. 
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Figure 9: Measure G{I) of seismic clustering in our BDE model of colliding 
cascades; see Eq. (1261) . The three curves correspond to the three synthetic 
sequences shown in Fig. [71 Reproduced from Zaliapin et al. [7], with kind 
permission of Springer Science and Business Media. 
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Figure 10: Regime diagram in the (A^,, A//) plane of the loading and healing 
delays. Stars correspond to the sequences shown in Fig. [71 Reproduced from 
Zaliapin et al. [7], with kind permission of Springer Science and Business 
Media. 
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Figure 11: Bifurcation diagram for the BDE seismic model: (a) rectangular 
path in the delay plane (A^, A^^); and (b) the measures Q and p, calculated 
along the rectangular path shown in panel (a). The transition between points 
(A) and (B), i.e. between regimes H and L, is very sharp with respect to 
the change in irregularity Q of energy release, while almost negligible with 
respect to the change in failure density p. The colored circles, triangles, and 
squares in panel (b) correspond to synthetic catalogs from regimes H, I, and 
L, respectively; these catalogs are produced for the points indicated along the 
rectangular path in panel (a), as well as for many scatter points that lie on 
a uniform grid covering the entire regime diagram, with the same resolution 
in Sh and 6l as those along the path. 
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Figure 12: Synthetic sequences corresponding to the key points along the 
rectangular path in parameter space of Fig. [TTh . The panels illustrate the 
transitions between the regimes H and L — panels (A) and (B); L and I — 
(C) and (D); and I and H — (E) and (F). The transition from (A) to (B) is 
very pronounced, while the other two transitions are smoother. Reproduced 
from Zaliapin et al. [7] , with kind permission of Springer Science and Business 
Media. 
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Figure 13: Solutions of the "partial BDE" flMl) : (a) for a single nonzero site 
at t = 0; and (b) the collision of two "waves," each originating from such a 
site. For the space and time steps 6z = Of = 1, this BDE is equivalent to the 
elementary cellular automaton (EGA) with rule 150; empty sites {ui{j) = 0) 
in white and occupied sites {ui{j) = 1) in black. 
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Figure 14: The solution of the BDE flM|) starting from a random initial state 
of length = 100. The qualitative behavior is characterized by "triangles" 
of empty {ui{i) = 0) or occupied = 1) sites but without any recurrent 

pattern; this behavior does not depend on the particular random initial state. 
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